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ABSTRACT 

We have developed a novel computer code designed to follow the evolution of cosmic-ray modified shocks, 
including the full momentum dependence of the particles for a realistic diffusion coefficient model. In this form 
the problem is technically very difficult, because one needs to cover a wide range of diffusive scales, beginning 
with those slightly larger than the physical shock thickness. With most finite difference schemes for Euler's 
equations the numerical shock thickness is at least one zone across, so this provides a lower bound on the physical 
scale for diffusive transport computation. Our code uses sub-zone shock tracking (LeVeque and Shyue 1995) and 
multi-level adaptive mesh refinement (Berger and LeVeque 1998) to provide enhanced spatial resolution around 
shocks at modest cost compared to the coarse grid and vastly improved cost effectiveness compared to a uniform, 
highly refined grid. We present and discuss the implications from our initial results. 

Subject headings: galaxy: globular clusters: general - hydrodynamics - ISM: supernovae remnants 



1. INTRODUCTION 

Diffusive shock acceleration (DSA) is now widely accepted 
as the model to explain the production of cosmic rays (CR) in a 
wide range of astrophysical environments (Drury 1983; Bland- 
ford and Eichler 1987; Berezhko and Krymskii 1988). The 
concept behind DSA, first-order Fermi acceleration of charged 
particles trapped between convergent flows across a shock, is 
quite simple. However, the full DSA problem is actually ex- 
tremely complex, because the nonlinear interactions between 
energetic particles, resonantly scattering waves and the under- 
lying plasma can become dominant effects. Important conse- 
quences of nonlinear interactions include such things as gen- 
eration and damping of the scattering wave field, injection of 
suprathermal particles into the CR population, as well as heat- 
ing and compression of the plasma flow due to the CR pres- 
sure. Owing to these complex nonlinear physics involved in the 
model, numerical simulations have been quite useful and suc- 
cessful in understanding the details of the acceleration process 
and dynamical feedback of the CRs to the underlying plasma 
(Falle and Giddings 1987; ElHson et al. 1990; Dorfi 1990; Kang 
and Jones 1991; Berezhko et al. 1994; Berezhko and Volk 
2000). 

In continuum approaches to numerical simulations of DSA 
theory, the CR diffusion-convection equation is solved at each 
of a large number of suprathermal momentum values simulta- 
neously with a set of fluid equations describing the flow as- 
sociated with the bulk, thermal plasma, including the nonlin- 
ear interactions between the plasma, CRs and scattering waves. 
Particle acceleration is effected by diffusion across velocity gra- 
dients in the motion of the scattering centers, which are usually 



assumed to be tied to the bulk flow. Pressure by the diffusing 
CRs, in turn, decelerates and compresses flow into the shock, 
forming a shock "precursor". Since that development elimi- 
nates the original, simple velocity jump seen by the CRs, the 
DSA is then modified according to details of the flow within the 
precursor, whose scales are characterized by the so-called dif- 
fusion length, D^iffip) = n{p) / u, where k is the spatial diffusion 
coefficient for CRs of momentum p, and u is the characteris- 
tic flow velocity against which the CRs must swim, e.g., (Kang 
and Jones 1991). Accurate solutions to the CR diffusion- 
convection equation require a computational grid spacing sig- 
nificantly smaller than D^iff, typically. Ax ^ 0.05Z)diff(/5). In 
a realistic diffusion transport model, it is thought that the dif- 
fusion coefficient should have a steep momentum dependence, 
k(p) oc p\ with s ^ 1-2. For the lowest energy CR particles 
the diffusion lengths (Z)diff(p)) are only slightly greater than the 
shock thickness, while they can be many orders of magnitude 
greater than that for the highest energy particles. Thus, a wide 
range of length scales is required to be resolved in order to solve 
the diffusion convection equation correctly for the model with 
a realistic diffusion coefficient. Previous numerical simulations 
which adopted the traditional flux-differencing method on a 
uniform grid were often forced to assume a weak momentum 
dependence, for example, s = 0.25 in Kang and Jones (1991). 

To overcome this numerical problem, Berezhko et al. (1994) 
introduced a "change of variables technique" in which the 
radial coordinate is transformed into a new variable, x(p) = 
exp[-(r-Rs)/Diiff(p)] where Rs is the shockradius, defined for 
each particle momentum for the upstream region. A uniform 
grid is used for the downstream region. This allowed them 
to solve the coupled system of gasdynamic equations and the 
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CR transport equation even when the diffusion coefficient has 
a strong momentum dependence {e.g., nip) ex. p). Their code is 
designed for simulations of supernova remnants, which are rep- 
resented by piston-driven spherical shocks in one-dimensional 
geometry. It is different from conventional Eulerian codes in 
several ways. Both gasdynamic equations and the CR trans- 
port equation are solved separately either side of the gas sub- 
shock. Then the gasdynamic solutions at both sides of the 
subshock are used to solve the Riemann problem, which de- 
termines how the subshock evolves. Also an iteration scheme 
is apphed to match the downstream and upstream solutions 
for the CR diffusion-convection equation at the subshock. In 
any case this has enabled them to explore several important is- 
sues regarding the particle acceleration at supernova remnants 
more fully than was possible before, e.g.,Berezhko et al. (1995, 
1996); Berezhko and Volk (2000). However, no consistency- 
check for this method has been attempted so far, since no exist- 
ing conventional codes can handle such a strongly momentum 
dependent diffusion coefficient. 

Fermi shock acceleration affects those particles with a mean 
free path greater than the shock thickness that can resonantly 
scatter with self-generated Alfven waves. In the so-called "ther- 
mal leakage" type injection model, the diffusion and accel- 
eration of these particles out of the suprathermal tail of the 
Maxwellian distribution determines the CR injection rate (Elli- 
son and Eichler 1984; Kang and Jones 1995). A self-consistent, 
analytic and nonUnear model for ion injection based on the 
interactions of the suprathermal particles with self-generated 
magneto-hydrodynamic waves in strong shocks has been pre- 
sented by Malkov (1998). By adopting this analytic solution, 
Gielseler et al. (2000) have developed a numerical treatment 
of the injection model at a strong quasi-parallel shock, which 
is then incorporated into the combined gas dynamics and the 
CR diffusion-convection code. Since the suprathermal parti- 
cles have mean free paths a few times that of thermal particles, 
resolving these smallest scales is of critical importance in esti- 
mating the injection and acceleration efficiency in such numer- 
ical simulations of the CR modified shocks. In fact, Gielseler 
et al. (2000) were able to run their simulations, with a con- 
ventional Eulerian scheme on a uniform grid, only up to the 
time when the maximum accelerated momentum was of or- 
der of Pmax/inpC ^ 1 for a Bohm type diffusion model because 
of severe requirements for computational resources needed to 
evolve the CR distribution to highly relativistic momenta. This 
calls for an alternative method comparable to Berezhko's code, 
which solves the CR diffusion-convection equation on a grid 
whose spacing scales with the diffusion length at each momen- 
tum value sampled. 

In this contribution, we present a new numerical scheme 
that follows CR modified shocks in one dimensional, plane- 
parallel geometry. We take advantage of the fact that the diffu- 
sion and acceleration of the low energy particles are important 
only close to the shock owing to their small diffusion lengths. 
They are simply advected along with the underlying gas flow 
far upstream and downstream of the shock. Thus it is neces- 
sary to resolve numerically the diffusion length of the particles 
only around the shock. So we first implement a shock track- 
ing scheme to locate the shock position exactly and then in- 
crease the grid resolution only around the shock by applying 
multi-levels of refined grids. Toward this end, we have adopted 
the shock tracking method of LeVeque and Shyue (1995) and 
the Adaptive Mesh Refinement (AMR) technique of Berger and 



LeVeque (1998), and modified the code to use multiple levels 
of grid refinement only around the shock. 

In the following section we outline our numerical methods, 
while in §3 we present and discuss our test results. Section 4 
provides a summary. 

2. NUMERICAL METHOD 

The diffusive transport model for CR acceleration separates 
the plasma into two components distinguished by scattering 
length. The bulk plasma consists of thermal particles whose 
scattering lengths are small enough to fit within a dissipa- 
tive shock. They are described by the standard gasdynamic 
equations with CR pressure terms added (McKenzie and Volk 
1982). The diffusion-convection equation, which describes 
the time evolution of the CR distribution function f{p,x,t) 
(e.^.,SkilUng (1975)) is given by 

^ = hv-ii)p^ + V-(Mx,p)Vf), (2-1) 
at 3 op 

where d/dt is the total time derivative in the fluid frame and 
the diffusion coefficient K(x,p) is assumed to be a scalar. As 
in our previous studies, the function g(p) = p'^fip) is solved in- 
stead of fip). Except for the special shock tracking and AMR 
features, our treatments of the underlying gas dynamics and 
the CR transport are relatively standard (Kang and Jones 1991; 
Gielseler et al. 2000), so we do not repeat them here. 

The spatial diffusion coefficient can be expressed in terms of 
a mean scattering length. A, as K{x,p) = | Av, where v is the par- 
ticle speed. The scattering length. A, and thus k(x,p), should be 
in principle determined by the intensity of resonantly interact- 
ing Alfven waves. For example, the Bohm diffusion model rep- 
resents a saturated wave spectrum and gives the minimum diffu- 
sion coefficient as kb = l/^^'gV when the particles scatter within 
one gyration radius (rg) due to completely random scatterings 
off the self-generated waves. This gives kb oc p^/ip^+ 1)'^^. 
Hereafter we will express particle momenta in units mpC. We 
consider here only the proton CR component. For our test runs, 
we will also adopt a power-law form as k(p) oc p" for low mo- 
menta (p < 1) in some models in addition to Kg. We note that 
the Bohm diffusion coefficient becomes K{p) cx p^ in the limit 
of p « 1 and k(p) oc p in the hmit of p » 1. In order to 
model ampUfication of self-generated turbulent waves due to 
compression of the perpendicular component of the magnetic 
field, the spatial dependence of the diffusion is modeled as 

k(x,p) = k(p)(pi/p(x)), (2-2) 
where pi is the upstream gas density. This form is also required 
to prevent the acoustic instability of the precursor (Drury and 
FaUe 1986; Kang, Jones and Ryu 1992). 

We also adopt the thermal leakage type injection model in- 
troduced in Kang and Jones (1995). In this model, below a 
certain momentum, pi, chosen high enough to include most of 
the postshock thermal population, the distribution is forced to 
maintain a Maxwellian form consistent with the local gas tem- 
perature and density determined from the gasdynanaical equa- 
tions. Above pi particles are allowed to evolve according to the 
diffusion-convection equation, so the form will deviate from 
Maxwellian. However, only for p > /72 > pi are they included 
in calculations of CR pressure and energy. We relate pi and 
P2 to the peak of the postshock Maxwellian distribution, /?,/„ as 
Pi = cipth, and p2 = C2Pth, and we assume ci = 2.5 and C2 = 3.0 
for all test runs here. Here pth corresponds to the peak in the 
partial pressure of thermal particles. The choice of p\ influ- 
ences the injection rate directly, since it determines the fraction 
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of suprathermal particles in the Maxwellian tail that can be in- 
jected into the CRs. 

2.1. Shock Tracking Method 

The hydrodynamic conservation equations are solved in the 
ID plane-parallel geometry by the wave-propagation algorithm 
described in LeVeque (1997). In this method a nonlinear Rie- 
mann problem is solved at each interface between grid cells, 
and the wave solutions (i.e. , speeds of waves and jumps asso- 
ciated with three wave modes) are used directly to update the 
dynamic variables at each cell. Within this method a sub-zone 
shock-tracking algorithm of LeVeque and Shyue (1995) can be 
incorporated easily, since the Riemann solutions tell us exactly 
how the waves propagate. The underlying Eulerian grid, which 
is called the "base" grid through this paper, has uniform cells. 
An additional cell boundary is introduced at the location of the 
shock, subdividing a uniform cell into two sub-cells. In the next 
time step, this cell boundary (shock front) is moved to a new lo- 
cation using the Riemann solutions at the current shock location 
(i.e. , x^^^ =x^ + Vs Af ) and the waves are propagated onto the 
new set of grid zones. Since the new grid is chosen so that the 
shock wave coincides exactly with an irregular cell boundary, 
the shock remains as an exact discontinuity without smearing. 
One advantage of using the wave-propagation method for the 
shock tracking scheme is that the large time step satisfying the 
Courant condition for the uniform grid can be used even if the 
shock is very close to the boundary of the uniform cell and so 
the sub-cell is much smaller than the uniform cell. 

The CR diffusion-convection equation is solved in two steps: 
1) the diffusion term is solved by the Crank-Nicholson scheme 
as described in Kang and Jones (1991). 2) the advection tennis 
solved by the wave-propagation method as for the gasdynamic 
variables. 

2.2. Adaptive Mesh Refinement 

Ideal gasdynamic equations in ID planar geometry do not 
contain any intrinsic length scales to be resolved, but once the 
precursor due to the CR pressure modification becomes signif- 
icant, the grid spacing should be fine enough to resolve the 
precursor structure. According to previous numerical studies 
e.g., (Jones and Kang 1990; Kang and Jones 1991), conver- 
gence of numerical solutions to CR modified shocks, especially 
the subshock transition, requires that the precursor be resolved 
with sufficient accuracy. While the full extent of the precursor 
increases with the CR pressure and is related with the diffusion 
length of the maximum accelerated momentum pmax, the dom- 
inant scale length of the precursor is similar to an averaged dif- 
fusion length of the particle populations with the greatest con- 
tribution to the CR pressure. Typically a strong shock becomes 
significantly modified due to nonlinear feedback from the CR 
pressure when the maximum acceleration momentum becomes 
mildly relativistic {i.e. pmax ~ 1). Thus in order to follow the 
development of the precursor and the time evolution of the CR 
modified shock, the gasdynamic equations should be solved on 
a base underlying grid whose spacing is smaller than the diffu- 
sion length of mildly relativistic particles. As discussed earher, 
it would be most natural to solve the CR transport equation on 
a grid whose spacing scales with the particle's momentum as in 
Berezhko et al. (1994). In that case, the CR distribution fix, p) 
should be mapped onto the base hydrodynamic grid in order to 
calculate the CR pressure and its dynamical feedback on the 
dynamics of the underlying flow. 



2.2.1. Laying down multi-level grids 

Here we take a different approach in which the immediate 
upstream and downstream regions around the shock are refined 
by applying multi-level grids with increasingly finer resolution 
by an integer factor of two, so that the transport of low energy 
particles right adjacent to the shock is at the most refined grid. 
Here we refer each level grid by the grid level index Ig which 
runs from to /max, corresponding to the base grid and the finest 
grid, respectively. Since the grid spacing decreases by an inte- 
ger factor, we can lay down the refined grid in such a way that 
cell boundaries ahgn between two adjacent levels. This fea- 
ture allows us to use a much simpler mapping scheme between 
two adjacent levels, compared to the case where non-integer 
refinement factors are used. We adopt the Adaptive Mesh Re- 
finement technique developed by Berger and LeVeque (1998). 
In the general version of the AMR code of Berger and LeV- 
eque (1998), the code identifies the refinement regions where 
the desired level of numerical accuracy is not achieved and the 
multi-level grids are generated within the refinement regions. 
Compared to that general version, a much simpler scheme is 
sufficient for our needs, since we only need to refine the region 
around the shock whose location is exactly known in our shock- 
tracking code. In fact, it is crucial to have the shock-tracking in 
order to lay down the multi-level grids around the shock, so that 
the shock remains near the middle of the computational domain 
at all levels. 

A fixed number of cells around the shock (A^rf) are identified 
as the "refinement region" on the base grid {i.e. Ig = grid). 
The 1st level grid is generated by placing 2Nrf cells within the 
refinement region, so each cell is refined by a factor of two. We 
use the refinement factor of two, since it is relatively simpler in 
terms of programming and it improves robustness and stability 
of the code. Then Nrf cells around the shock out of 2Nrf cells 
on the 1st level grid are chosen to be refined further to the 2nd 
level grid, making the length of the refinement region a half of 
that in the 1st level grid. Here the refinement region is chosen 
so that the shock is always in its middle. The same refinement 
procedure is apphed to higher level grids. So at all levels, there 
are 2A^rf cells around the shock, but the length of the computa- 
tion domain is shrunk by a factor of 2 from the previous level. 
At each level the grid spacing is given by Ax{lg) = Ax(0)/2's, 
where Ax(0) is the grid spacing at the base grid. Fig. 1 shows 
an example of refined grid levels up to /max = 2 with A^ff = 4. 
The value of A'rf should be chosen so that the refined region at 
the base grid includes most of the precursor during the early 
stages when particles get accelerated to mildly relativistic ener- 
gies. In real test simulations presented here, typical values of 
Nrf = 100-200. The spatial extent of the highest refined grid 
can be much smaller than the length scale of the precursor, since 
it need only resolve structures "seen" by the lowest energy CRs 
immediately next to the shock. 

In order to ensure that the shock remains near the middle of 
the computational domain at all grid levels during the time in- 
tegration of one time step of the base grid, we do the following 
procedure. First, the velocity in the refined grid is transformed 
so that the shock is at rest in the frame of the numerical simula- 
tion at each level except the base grid. Secondly, the multi-level 
grids are redefined at each time step so that they center around 
the new location of the shock, and all hydrodynamic variables 
and the particle distribution function g{p) are mapped onto the 
newly defined grids. Thus the refinement region at all levels 
is moving along with the shock. Although the original shock 
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tracking code of LeVeque and Shyue (1995) can treat multiple 
shocks, we modify it in the current version of our code to in- 
clude only one shock, in order to keep the code structures as 
simple as possible for our initial studies. 

2.2.2. Time integration 

Integration of the gasdynamic variables and advection of 
the CR distribution function are done by the wave-propagation 
method at each level grid. The time step at each level, Af(Zg), 
is determined by a standard Courant condition, that is, Af (Z^) = 
Q3Ax{lg)/ maMu+Cs), where u and Cs are the flow velocity and 
sound speed at each cell, respectively. If the highest level is 
specified to be Z^ax = 1, to advance from f" to f"+' = f" + Af(0) 
at the base grid, we need the following steps: 1) all equations 
are integrated at the base grid. 2) the same is done twice with 
Af(l) = Af(0)/2 at the Ig = 1 level for the refined region. The 
cells immediately outside the refinement region at the Ig = Q 
level provide the necessary boundary conditions for the integra- 
tions. Here we need boundary conditions at two spatial points 
at both ends of the refined grid, qq, qi, qiN^f, and q2Nrf+i< and 
also at two points in time, t" and f""""'/^. Some of them are in- 
terpolated in time and space coordinates from the variables de- 
fined at the grid one level below. 3) the values at the /g = 1 
level are mapped onto the refinement region at the base grid. 
4) finally the values at the interface just outside the refinement 
region at the base grid should be corrected to preserve global 
conservation. The base idea and also the detailed procedure ap- 
plied for one level of refinement can be found in Berger and 
LeVeque (1998). When /max is greater than 1, the same proce- 
dures should be repeated recursively at each time step at each 
level. The variables within the refinement region at the Ig level 
grid are replace with the more accurate values at ihelg+l level 
grid after the corresponding pair of time steps at the /g -I- 1 level 
are completed. In one time step in the base grid, we carry out 
2's time steps at the Ig level grid, so total number of time steps 
in all level grids required to advance one time step in the base 
grid becomes the sum of 2 -i- 2^ -i- . . . -i- 2'""« . 

Once advection of the CR distribution is done, diffusive 
transport including the first-order acceleration is solved in a 
separate step. Although an implicit Crank-Nicholson scheme is 
used for the diffusion term, the time step is still restricted by the 
acceleration term as AtcN(lg) = 0.5min[3A(log/7)/(9M/9x),]. 
Within a single hydrodynamic time step, several Crank- 
Nicholson time steps are performed if At{lg) > Atcyiilg). Num- 
ber of this sub-cycling is about 4-5 for the momentum bin size 
considered here (i.e. , Ailogp) = 0.026). 

As we go up and down on the ladder of "time stepping" for 
one time step in the base grid, the base grid being the lowest, the 
values at a coarser grid propagate upward as boundary condi- 
tions for a finer grid and more accurate values at the finger grid 
propagate downward by being mapped onto the coarser grid. 

3. TEST RESULTS 

In this section we present some test simulations with our 
CR/AMR code. The dynamics of the CR modified shock de- 
pends on four parameters: the gas adiabatic index, jg = 5/3, 
gas Mach number of the shock, M = V^/cs, f3 = Vs/c, and the 
diffusion coefficient, where Cj and c are the upstream sound 
speed and the speed of light, respectively. For all simula- 
tions we present here, f3 = 10"^, and 7g = 5/3. We considered 
three values for Mach number, M = 5, 10, and 20 for the ini- 
tial shock jump by adjusting the preshock gas pressure. The 



initial jump conditions in the rest frame of the shock for all 
test problems are: pi = l.,Pg,i = 1.5 x 10~^(20/Mf,ui = -1. 
in upstream region and p2 = 4.,/'g,2 = 7.5 x 10"', M2 = -0.25 
in downstream region. Here the velocities are normalized to 
the initial shock speed, = 3000kms"'. Normalization of the 
length and the time variables depends on the diffusion coeffi- 
cient: Kip) = Kphys/Ko, where Kip) is the computational coeffi- 
cient, Kphys is the physical value, and Kg is the normalization 
constant. So the corresponding normalization constants are: 
to = Ko/V^, rg = Ko/Vs. These correspond to roughly the dif- 
fusion length and diffusion time scale for p ~ 1 . For the Bohm 
diffusion coefficient with the magnetic field of 1 microgauss, 
for example, kb = p^/ip^+ 1)'^^ with Hg = 3.13 x lO^^cm^ s"\ 
so that to = 3.5 X lO^s and rg = 1.05 x lO'^^cm. The particle 
number density, rig is arbitrary, but the gas density and pressure 
are normaUzed to pg = mpHg and Pg = PoV^, respectively. 

We use 230 logarithmic momentum zones in log(p)=[- 
3.0,+3.0] and the momentum is in units of iripC. The distri- 
bution function fip) is expressed in units of /„ = ng/im^c)^, so 
that 4-n- J f p^ dp = p. 

3.1. Test of Refinement 

In this section we consider the M = 20 shock with a diffu- 
sion model Kip) = p. In order to see how the CR/AMR code 
performs at different resolutions, we ran the simulations with 
different levels of refinements, /max = 0, 1,2,3,4. The numeri- 
cal domain is [-25,+751 and the number of cells, = 2000, so 
the grid spacing is Ax(0) = 0.05 which corresponds to 1/20 of 
the diffusion length of the particles of p = I. The number of 
refined ceUs around the shock is Nr/ = 200 on the base grid and 
so there are 2Nrf = 400 cells on each refined level. Since the 
diffusion length of pi ^ 0.01 is Ddiff = 0.01, the transport of the 
suprathermal particles can be resolved at the /max = 4 level. 

With the refinement of /max = 4 level, the CR transport for 
the particles at the injection pool should be resolved and so the 
evolution of the CR should be converged. Fig. 2 shows the time 
evolution of the model shock with the maximum refinement of 
/max = 4. This shows how the precursor develops and modifies 
the shock structure as the CR pressure increases in time. The 
numerical frame is chosen so that the initial shock moves to the 
right with Us = 0.05, but the simulated shock drifts to the left 
due to the CR pressure . Fig. 3 shows the density structure at 
each refined grid (soUd lines) for /g = 1 , 2, 3 , and 4 levels, super- 
posed on the density structure at the base grid (dotted lines) in 
the /max = 4 simulation shown in Fig. 2. This demonstrates how 
the size of refinement region decreases at higher levels and how 
the refinement regions move along with the shock. As indicated 
by two points at downstream and upstream of the subshock in 
the top-left panel, the shock is tracked as a perfect discontinuity. 

Given the same resolution at the base grid, the simulations 
with larger refinement levels show faster acceleration and faster 
growth of Pc. Fig. 4 shows how the particle distribution 
(Eip) = fip) P'^) at the shock evolves with time in the simula- 
tions with different refinement levels for the same shock model 
shown in Fig. 2: there are 5 curves corresponding to /^ax = 0, 
1, 2, 3, and 4 in each panel. They show the typical Maxwelhan 
distribution that peaks at pih ~ 10"^^ and the CR distribution 
that asymptotes to a power-law as time increases. For /^ax = 
(dotted lines), the cell size Ax = 0.05 is too large for the diffu- 
sion of the particles in the injection pool to be treated correctly, 
so the injection and the acceleration are under-estimated and 
the slope just above the injection pool is steeper than the canon- 
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ical strong-shock test-particle spectrum of f{p) ^ p'"^. For the 
highest refinement case, /max = 4 (solid lines), the CR pressure 
becomes dominant over the gas pressure and the compression 
ratio at the subshock becomes 3.3. As a result, the distribution 
function steepens from a power-law of p~'^ at lower momentum 
just above the injection momentum. But it flattens at higher 
momenta since high energy particles diffuse on a much larger 
scale and sample a larger velocity jump. These opposite trends 
lead to a concave curve in the middle, which is a typical sig- 
nature of nonlinear effects (Berezhko and Ellison 1999). Fig. 
5 shows how the CR pressure increases with time in the sim- 
ulations with different refinement levels shown in Fig. 4. It 
demonstrates that injection and acceleration are much slower 
in under-resolved simulations (/max < 4). From comparison be- 
tween results of /^ax = 3 (dashed line) and /^ax = 4 (solid lines), 
we have concluded that the simulation is mostly converged for 
/max = 4. For example, the CR pressure at the shock at ^ = 50 and 
60 is the same for Z^ax = 3 and 4 simulations, although the shock 
position is slightly different in the two cases due to slightly dif- 
ferent evolution in early stages. 

The required computing time increases with /^ax and Nr/, 
given the same resolution at the base grid. For the simulation 
considered here in which Nrj/N = 0.1, the computing time in- 
creases by factors of 1.5, 2.3, 4, 7 for the maximum refinement 
levels /max = 1, 2, 3, 4, respectively, compared with the case of 
no refinement (Z^ax = 0). The computing time would increase 
by factors of (2'"")^, if the simulations were done on a uni- 
form grid spacing that matches the cell size at the highest re- 
fined level grid. Fig. 6 compares the computing time for these 
two cases. Since only the precursor region needs to be refined, 
our CR/AMR code will be most cost-effective for simulations 
where the precursor is only a small fraction of the total compu- 
tational domain. 

3.2. Convergence Test 

In this section we explore how a simulation with our 
CR/AMR code using multi-level refinements would be com- 
pared with that with no refinement but on a single uniform grid 
of the same spacing as the highest level of the other. For this 
test, we use the following diffusion model: K,\=p/ \/2 for p <l 
and Ath = p^/(p^ + 1)^''^ for p>l. This model has a Bohm- 
type diffusion at higher momenta, but much weaker momen- 
tum dependence at lower momenta. This allows us to use larger 
grid spacing near the injection pool {p ~ 0.01) compared with a 
Bohm diffusion model which scales as p^ for p «\. The ini- 
tial shock parameters are the same as the M = 20 shock model 
in previous section. We consider the following three models: 
Model A, a uniform grid with Ax(0) = 0.1 and /max = 0, Model 
B, an adaptively refined grid with Ax(0) = 0.1 and /max = 5, and 
Model C, a uniform grid with Ax(0) = 0.1/32 and /max = 0. 
The numerical domain is [-50,H-50] and the number of cells, 
N = 1000 for Models A and B, and = 32000 for Model C. So 
Model C has the same grid spacing as the /max = 5 level grid of 
Model B. The computing time for Model C with a single fine 
grid is about 150 times longer than that for Model B with the 
refinement. We expect that in the Model A simulation the trans- 
port of low energy particles are under-resolved, so it is included 
only for comparison. We choose N,-/ =100 for Model B, so 10 
% of the base grid is refined. The major difference between 
Models B and C is that the refined grid at /max = 5 covers only 
about 1/32 of the precursor in Model B, while the entire grid of 
Model C has the finest resolution. 



Fig. 7 shows the evolution of the gas density and the CR pres- 
sure, and the CR distribution function at the shock in Models B 
and C. In Model B with refinement the CRs are accelerated a 
little faster than in Model C with a single fine grid, so that the 
difference in the CR pressure at the shocks is 8.3 % f = 20, but 
this fraction decreases to 4.6 % at f = 40 and to 3.0 % at f = 60. 
Thus, the two methods appear to converge at slightly different 
rates, but both give reasonable results once the resolution next 
to the shock is refined to resolve £>diff(/'i) by more than an or- 
der of magnitude. This test convinces us that our CR/AMR 
code can perform the intended simulations with a reasonable 
accuracy in a very cost-effective way. 

3.3. Dependence on Diffusion Model 

Finally in this section we explore briefly how different dif- 
fusion models affect evolution of the injection and the acceler- 
ation efficiency in CR modified shocks. First we consider the 
M = 20 shock and the following three diffusion models. For 
lower momenta, p < 1, Model Kl uses k\ i = p/V2, Model 
K2 Ki,2 = p' VV^, and Model K3 K13 = pV(p2 + l)i/2. All 
three models use the Bohm type diffusion at higher momenta, 
that is, Kh = p^ /{p^ + 1)'/^ for p > I and so k\ continuously 
matches onto Kh at p = 1. The grid spacing for the base grid, 
Ax(0) = 0.05,0.025, and 0.005 for Models Kl, K2, and K3, re- 
spectively. The maximum refinement level, /max = 4,7, and 7 
for Models Kl, K2, and K3, respectively. These parameters are 
chosen so that the grid spacing at the finest grid is fine enough 
to treat the low momentum particles near pi with the assumed 
diffusion model. While Model Kl was integrated for t = 100 
and Model K2 for t = 60, Model K3 is integrated only up to 
t = 20 due to much longer required computing time. 

Fig. 8 shows the gas density, the CR pressure, and the CR 
distribution function at the shock for the M = 20 shock at f = 20 
simulated with the three different nip) (Kl: soUd line, K2: dot- 
ted line, K3: dashed line). Model K3 has the smallest D^{p), 
while Model Kl has the largest D^iffip) for the particles in the 
injection pool (p ^ 0.01). So the injection takes place the fastest 
and the CR pressure increases most efficiently in Model K3 
during the early evolution. According to earUer evolution (not 
shown), the time evolution of these three models differs sig- 
nificantly when the maximum accelerated momentum is still 
nonrelativistic (pmax << 1)- At f = 20, however, pmax ^ 2 and 
mildly relativistic particles dominate the CR pressure, so the 
different diffusion models at nonrelativistic momenta no longer 
play a significant role. At this time all three models evolve in 
a similar way, since the same Kh is used for all models when 
p> I. Especially Models K2 and K3 show very similar evolu- 
tion up to this time. Subsequent development of the shocks is 
almost independent of the form for ki. The bottom right panel 
shows the CR distribution at the shock at f = 60 in Models Kl 
and K2 (K3 was ended at f = 20), demonstrating the very similar 
CR distribution evolution in Models Kl and K2 at later times. 
Considering the earlier trend that Models K2 and K3 are al- 
ready very similar at t = 20, we can deduce all three models 
evolve the same way when /^max becomes much larger than one. 
This imphes that k (x p can be used instead of kb as long as 
the detailed evolution at early stage when the particles are still 
mostly nonrelativistic is not taken seriously. Simulations with 
k(p) oc p model allow much coarser grid spacings to follow the 
lowest momentum particles than those with a Bohm type dif- 
fusion, which reduces the required level of refinements and the 
associated costs. The same set of simulations were repeated for 
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M =10 and M = 5 shocks, and we came to the same conclusion 
(see Fig. 9 for results from the Mach 10 simulation). This val- 
idates the notion that K{p) oc p can be used instead of Kb in the 
CR/hydro simulation of SNRs (Berezhko et al. 1995; Berezhko 
and Volk 2000). 

We also ran a model with a pure power-law type diffusion 
model, that is, k = p for all momenta until ? = 10'' to study 
the long-term evolution. For this simulation, the computational 
domain is increased to [-4000,h-4000], the zone number in the 
base grid to N = 40000 and the maximum refinement level to 
^max = 6. Fig. 10 shows the shock structure, the CR distribu- 
tion at the shock, and its power-law slope, q = -{d\af /d\n p), 
at t = 200, 800, 2 x 10^ 5 x 10^ and lO'^. The CR pressure 
at the shock seems to have reached a quasi-steady state value 
at t > 200, even though the maximum momentum continues 
to increase with time. The compression ratio at the subshock 
is ~ 3.1 at ? = 10^ which leads to the test-particle slope of 
q ^ 4.5. So the particle distribution near the injection pool re- 
flects this slope. Although the CR pressure at the shock stays 
more or less constant after t = 200, and the total shock jump 
does also, the precursor and its associated velocity structure 
broadens with time. Thus, since Kip) does not change over time 
in this idealized simulation, the particles of a given momentum 
sample a smaller velocity jump, Auip), around the shock, as 
the precursor broadens, resulting in slightly steeper slope. On 
the other hand, the highest momentum particles sample some- 
thing close to the full velocity jump on their diffusion scales, 
so the slope flattens gradually to q ^ 3.2 toward pmax, which 
corresponds to the total compression ratio, rtot ~ 11-12. This 
hardening of the CR distribution from ^ ~ 4.5 to ^ ^ 3.2 pro- 
duces concave curves in log(^ = fp'^) versus log(p) plot. This 
illustrates the importance of following correctly the non-linear 
feedbacks between the CRs and the dynamics inside the pre- 
cursor, which requires one to resolve numerically all relevant 
scales. 

Finally, Fig. 1 1 shows how the postshock gas temperature 
decreases as the CR pressure becomes dominant in the precur- 
sor and how the injection parameter settles down to a constant 
value at ry ~ 6 X 10~^ after the shock has reached a quasi- steady 
state. The injection parameter r] is defined in Kang and Jones 
(1995) and represents the fraction of the incident proton flux 
that is injected into the CR population at the shock. We also 
plot the adiabatic index of the CR population at the shock which 
also settles down to a constant value at 7c ^ 1.37. The injec- 
tion parameter and the CR adiabatic index along with the mean 
diffusion coefficient are three free parameters for the so-called 
two-fluid model for the CR modified shock simulations. 

4. SUMMARY 

We have developed a new hydro/CR dynamics code which 
is specifically designed to solve the time dependent evolution 
of CR shocks. Diffusive shock acceleration of the CR parti- 
cles depends on the diffusion of particles whose momenta span 
many orders of magnitude. Since the length and time scales for 
evolution of the CR kinetic equation scale directly with the dif- 
fusion coefficient, an accurate solution to the problem requires 
that one include all of those scales in the simulation, beginning 
just outside the gas subshock thickness. Thus, in order to fol- 
low accurately the evolution of a CR modified shock, it is nec- 
essary to resolve the precursor structure upstream of the sub- 
shock and, at the same time, to solve correctly the diffusion of 
the low energy particles near the injection pool. These low en- 



ergy particles have diffusion lengths that are much smaller than 
the scale length of the precursor, so a large dynamic range of 
resolved scales is required for CR shock simulations. To solve 
this problem generally we have successfully combined a power- 
ful "Adaptive Mesh Refinement" (AMR) technique (Berger and 
Colella 1989; Berger and LeVeque 1998) and a "shock track- 
ing" technique (LeVeque and Shyue 1995), and implemented 
them into a hydro/CR code based on the wave-propagation 
method (LeVeque 1997). The AMR technique allows us to 
"zoom in" inside the precursor structure with a hierarchy of 
small, refined grid levels applied around the shock. The shock 
tracking technique tracks hydrodynamical shocks and main- 
tains them as true discontinuities, thus allowing us to refine the 
region around the shock at an arbitrary level. The result is an 
enormous savings in both computational time and data storage 
over what would be required to solve the problem using more 
traditional methods on a single fine grid. 

The code has been applied to simulations of CR modified 
shocks with several diffusion coefficient models with strong 
momentum dependence, which were not possible previously 
due to severe computational requirements. The main conclu- 
sions from the simulations can be summarized as follows: 

1 . Our CR/AMR technique code proves to be very cost ef- 
fective. In typical simulations where 10% of the base 
grid is refined with /niax levels, for example, the comput- 
ing time increases by factors of (2'"'")"^ compared with 
the case of no refinement (/max = 0). It should be com- 
pared with the time increases by factors of (2'™")^ for the 
simulations of an uniform grid spacing that matches the 
cell size at the /^ax - th refined level grid. In a simula- 
tion where the precursor scale is only a small fraction of 
the computational domain, the advantage in computing 
time of the refined multi-level grid over a single fine grid 
becomes even greater. 

2. A convergence test is performed for a Mach 20 gas 
shock with Vg/c = 0.01, which evolves into a CR pres- 
sure dominated shock. Comparison between a simula- 
tion on a coarse grid with multi-level refinement and an- 
other simulation on a single fine grid without refinement 
has demonstrated that our CR/AMR code can perform 
the intended simulation with reasonable accuracy at a 
much shorter computing time. The difference in the CR 
pressure in two test simulations is around 10% in early 
evolution, but deceases to a few % after the shock has 
reached a quasi-steady state in later evolution. The re- 
quired computing time is reduced by a factor of 150 in 
the AMR simulation. 

3. We also carried out a set of simulations when three 
different diffusion models, ki = p/V2, /?'^/\/2, and 

//(/?2+ 1)1/2 foj. ^ ^ 

are included, while a Bohm 
type, Kh = p^/ip^+ 1)^''^ is assumed for p > 1. Three 
simulations generate similar results once the CR pres- 
sure is dominated by the relativistic particles (p > 1), 
when the maximum acceleration momentum becomes 
Pmax >> 1. Thus a diffusion model of k oc p can be 
used instead of a Bohm model as long as one does not 
focus on the early evolution when p^^^x < 1- Since we 
can use much larger grid spacings in simulations with 
K oc p model than those with a Bohm type diffusion at 
low momenta p « 1, the required level of refinements 
and so the computing resources can be reduced greatly. 
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4. For a Mach 20 shock, with an injection rate of 77 ~ 6 x 
lO""*, the shock becomes CR dominated and develops a 
significant precursor. Since the flow is decelerated grad- 
ually through the precursor, the velocity jump that the 
CR particles sample across the shock depends on the dif- 
fusion length of the particle, that is, Au(p) = fcti(K(p)). 
So the slope of the particle distribution function, defined 
as q{p) = -{d\nf/d]np), increases with p. In the simu- 
lated shock, the compression ratios across the subshock 
and across the total transition are 3.1 and 11, respec- 
tively, so fip) is (X p""*^ at low energy momenta but flat- 
tens to fip) oc p~^-^ at high energy momenta just below 
Pmax- This demonstrates that nonUnear feedbacks be- 
tween the precursor dynamics and the CR injection and 
acceleration should be treated accurately in numerical 
simulations of CR shocks. 



We are currently implementing a numerical scheme for the 
self-consistent injection model by Gielseler et al. (2000) which 
is based on the plasma-physics study of the nonlinear wave- 
particle interactions in shocks presented by Malkov (1998). 
This will allow us to eliminate any free parameters for the in- 
jection process from the CR shock simulations. We intend also 
to extend the code to treat spherical shocks in order to study CR 
acceleration in supernova remnants. 
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Fig. 1 . — Layout of the base grid and two refined grids. Here N^f = 4 cells around the shock are refined by a factor of two. The shock is indicated by the dotted 
lines. 
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Fig. 2. — Time evolution of the M = 20 sliock with /max = 4 refined grid levels at ( = 10, 20, 30, 40, 50, and 60. The shock moves to the left, so the right most 
plots correspond to the earliest time f = 10. The solid lines are for the refinement region at = 1 grid, while the dotted lines represent the shock structure on the base 
grid. 
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Fig. 4. — Distribution function (g = fp^) at the shock for the M = 20 shock at f = 10, 30, 50, and 60. The results of the simulations with the maximum refined 
grid level /max = (dotted lines), 1 (long dashed), 2 (dot-dashed), 3 (dashed) and 4 (solid) are plotted for comparison. 
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Fig. 5. — Same as Fig. 4 except CR pressure is plotted. 
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Fig. 6. — The solid line shows the ratio of computational time required to include Lmax levels relative to the time required vi'ith no refinement. Here the number 
of cells at the base grid is W = 2000 and N,.f = 200 cells around the shock are refined. The dotted line shows the same ratio for the case when the finest resolution is 
appUed over the entire grid, that is, (2'-">")^. 
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Fig. 7. — Time evolution of the shock structure and the distribution function g = fp^ at the shock for Models B (dashed line) and C (soUd Une). The shock moves 
to the left, so the right most plots correspond to the earliest time f = 20 in the gas density and CR pressure plots. 



Cosmic Ray Modified Shocks 



15 




-15 -10 -5 5 10 -15 -10 -5 5 10 




-3-2-1012 -3-2-1012 
log(p/mc) log{p/mc) 



Fig. 8. — Comparison of the M = 20 shock structure at J = 20 for Models Kl (soUd), K2 (dotted), and K3 (dashed) with different diffusion coefficients (top 
panels). Comparison of the distribution function g = fp^ at the shock at r = 20 for Models Kl, K2, and K3 (bottom left) and at r = 60 for Models Kl and K2 (bottom 
right). 




Fig. 9. — Same as Fig. 8 except that the Mach number is 10. 
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Fig. 1 0. — Time evolution of the flow velocity and the CR pressure for the M = 20 shock with K(p) oc p for all momenta. Also the CR distribution function g at 
the shock and its power slope g = -dhif/dlnp. The shock moves to the left, so the right most plots correspond to the earliest time t = 200 in the flow velocity and 
CR pressure plots. 
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Fig. 11. — Time evolution of the postshock shock temperature, injection rate, and the adiabatic index of the CRs at the shock from the same simulation shown in 
Fig. 10. 



